* This file creates the seasonal hunger tables in Figure 2.

foreach pop in "BGD" "NPL" {

**********

* Set colors and other graph attributes
graph set window fontface "Helvetica"
	local lscol `""230 159 0%15""'
	local hvcol `""0 158 115%15""'
	local cvcol `""0 114 178""'

if "`pop'"=="BGD" {
tempfile temp

* 1) Combine midline and endline hunger data

* Get the hunger (binary) data by month
use "$da/BGD_endline_clean.dta", clear
keep hhid village end_seasonality_*
rename end_* mid_*

append using "$da/BGD_midline_clean.dta", keep(hhid mid_seasonality_*) gen(midline)
rename mid_seasonality_* pre_hunger_*

* Reshape to long over BGD months (dec. 15--jan 15 = 9)
egen idno = group(hhid)
rename *_jaishthha *2
rename *_baishakh *1
rename *_chaitra *12
rename *_falgum *11
rename *_magh *10
rename *_poush *9
rename *_agrahayan *8
rename *_kartik *7
rename *_ashwin *6
rename *_bhadra *5
rename *_sraban *4
rename *_ashar *3

reshape long pre_hunger, i(idno midline) j(month)

* Collapse duplicate months asked in mid and endline 
gen byte series = (midline==0) | inrange(month,3,9)
// Trust the endline
duplicates tag hhid series month, gen(dup)
drop if midline & dup
drop midline dup

* Convert to USA months
replace month = mod(month-10,12)+1

* 2) Add in COVID-19 data
rename hhid caseid
merge m:1 caseid month using "$da/BGD_covid_clean.dta", assert(match master)
drop _merge

rename caseid hhid
rename weights cvd_weight
rename hungry cvd_hunger

* Clear out 2019 months
replace cvd_weight = . if series==0
replace cvd_hunger = . if series==0

* 4) xtset and make it pretty
egen hhid_no = group(hhid)
order hhid hhid_no village month series *_hunger cvd_weight

* 5) Run regressions (weighted by sampling likelihod for post-COVID)
quietly areg pre_hunger i.month, a(hhid_no) vce(cluster hhid_no)
predict b_pre_hunger if !missing(pre_hunger), xb
predict se_pre_hunger if !missing(pre_hunger), stdp
	
quietly areg cvd_hunger i.month if !missing(cvd_weight) [aw=cvd_weight], a(hhid_no) vce(cluster hhid_no)
predict b_cvd_hunger if !missing(cvd_hunger), xb
predict se_cvd_hunger if !missing(cvd_hunger), stdp
	
collapse b_* se_*, by(month) fast
sum b_*, d

* Renumber months to put June at start and end, bump by 0.5
replace month = mod(month-6,12)+0.5
expand 2 if month==0.5, gen(newmon)
replace month = month+12 if newmon
drop newmon 
sort month
gen m2 = month - 0.5 if inrange(month,7,12)

* Generate standard error bars and convert to percents
foreach vbl in pre_hunger cvd_hunger {
	replace b_`vbl' = 100*b_`vbl'
	gen hi_`vbl' = b_`vbl' + 196*se_`vbl'
	gen lo_`vbl' = b_`vbl' - 196*se_`vbl'
}

* Locals for graphing
	local top 25.7
	local bottom -2
	local Tht 24.7
	local ls1 3
	local ls2 5.25
	local lsgap .03
	local lst = (`ls1'+`ls2')/2
	local lsx1 = `ls1'+`lsgap'
	local lsx2 = `ls2'-`lsgap'
	local hv1 5.5
	local hv2 6.5
	local hvt = (`hv1'+`hv2')/2
	local cov 9
	
	local shading (scatteri `bottom' 5.5 `bottom' 6.5, recast(area) color(`hvcol') lwidth(none)) ///
			(scatteri `top' 5.5 `top' 6.5, recast(area) color(`hvcol') lwidth(none)) ///
			(scatteri `bottom' 9 `Tht' 9, recast(line) lcolor(maroon) lpattern(dot))
	local shade2 xline(3.03(0.03)5.21 , lcolor(`lscol') lwidth(thick) lpattern(dot)) ///
			xline(3(0.06)5.25 , lcolor(`lscol') lwidth(vthin) )
	local words text(`Tht' 4.125 "Lean" "Season") text(25.2 9.05 "COVID-19") text(`Tht' 6 "Rice" "Harvest")
	local arrows (pcarrowi 20 11.5 18.5 10.7 (2) "2020", color(`cvcol') mlabcolor(`cvcol') mlabsize(medium)) /// 
			(pcarrowi 8.8 11.4 4.5 11 (12) "Previous Year", color(gs0) mlabcolor(gs0) mlabsize(medsmall))

	local yaxis ytitle("Percent of Households that are Food Insecure", axis(1)) ylabel(#5,glcolor(gs12%25))
	local xaxis xlabel(1 "Jul" 2 "Aug" 3 "Sep" 4 "Oct" 5 "Nov" 6 "Dec" 7 "Jan" 8 "Feb" 9 "Mar" 10 "Apr" 11 "May" 12 "Jun", angle(45)) xtitle("")
}

if "`pop'"=="NPL" {
use "$da/NPL_panel_clean.dta", clear
append using "$da/NPL_recall_clean.dta", gen(recall)

egen rec_hh = max(recall), by(hhid)
keep if rec_hh == 1

* Rounds 1-7: 2019/2020 survey data; Rounds 9-20: recall data
_pctile date if round == 6, nquantiles(2)
replace round = 7 if date > `r(r1)' & !mi(date)
replace round = 8 + month if recall == 1

* Generate numeric code for hhid
egen hhno = group(hhid)
sort hhno
xtset hhno

***** Make data set for regression analysis
keep ___foodInsSub ___foodIns_r round recall date hhno

* Main effect for pre-data
xtreg ___foodIns_r i.round if recall==1, fe vce(cluster hhno)
predict b_pre_hunger if recall==1 & !missing(___foodIns_r), xb
predict se_pre_hunger if recall==1 & !missing(___foodIns_r), stdp

* Main effect for post-data
xtreg ___foodInsSub i.round if recall==0, fe vce(cluster hhno)
predict b_cvd_hunger if recall==0 & !missing(___foodInsSub), xb
predict se_cvd_hunger if recall==0 & !missing(___foodInsSub), stdp

collapse (median) date (mean) b_* se_*, by(round) fast
replace date = td(01jun2020) if round==7		// Looks better graphically
sort round

* Normalize each index so it has mean 0 and sd 1 in Sept-Decs
foreach s in pre cvd {
	sum b_`s'_hunger if date < date("Jan 1, 2020", "MDY") & date > date("Aug 1 2019", "MDY")
	local mn_`s' = `r(mean)'
	local sd_`s' = `r(sd)'
	replace b_`s'_hunger = (b_`s'_hunger - `r(mean)')/`r(sd)'
	replace se_`s'_hunger = (se_`s'_hunger)/`r(sd)'
}

*Bump back the "July 2020" recall observation to July 2019
replace date = date - 365 if round>7 & date > date("30 June, 2020", "DMY")

*Duplicate last recall observation and stick it at the beginning
sort date
quietly count
expand 2 in `r(N)', gen(newvar)
replace date = date - 365 if newvar==1
drop newvar
sort date
rename date month

* Generate standard error bars 
foreach s in pre cvd {
	gen hi_`s'_hunger = b_`s'_hunger + 1.96*se_`s'_hunger
	gen lo_`s'_hunger = b_`s'_hunger - 1.96*se_`s'_hunger
}

* Locals for graphing
	local top 2
	local bottom -1.24
	local Tht 1.86
	local ls1 = date("2 Aug, 2019", "DMY")
	local ls2 = date("1 Oct 2019", "DMY")
	local lst = (`ls1'+`ls2')/2
	local lsgap = 1
	local lsgap2 = 2*`lsgap'
	local lsx1 = `ls1'+`lsgap'
	local lsx2 = `ls2'-`lsgap'
	local hv1 = date("11 Oct 2019", "DMY")
	local hv2 = date("15 Nov 2019", "DMY")
	local hvt = (`hv1'+`hv2')/2
	local cov = date("5 Mar 2020", "DMY")
	local cvt = (`top'+`Tht')/2
	
	
	loc ax2 = month[17] + 5
	loc ay2 = b_cvd_hunger[17] + .05
	loc ax1 = `ax2' + 15
	loc ay1 = `ay2' + .20
	loc bx2 = month[16] + 5
	loc by2 = b_pre_hunger[16] + .05
	loc bx1 = `bx2' + 15
	loc by1 = `by2' + .20
	loc cx2 = month[12] + 5
	loc cy2 = b_cvd_hunger[12] - .05
	loc cx1 = `cx2' + 15
	loc cy1 = `cy2' - .20
	local start = date("2 July, 2019", "DMY")
	local end = date("1 June, 2020", "DMY")
	
	local shading (scatteri `bottom' `hv1' `bottom' `hv2', recast(area) color(`hvcol') lwidth(none)) ///
			(scatteri `top' `hv1' `top' `hv2', recast(area) color(`hvcol') lwidth(none)) ///
			(scatteri `bottom' `cov' `Tht' `cov', recast(line) lcolor(maroon) lpattern(dot))
	local shade2 xline(`lsx1'(`lsgap')`lsx2' , lcolor(`lscol') lwidth(thick) lpattern(dot)) ///
			xline(`ls1'(`lsgap2')`ls2' , lcolor(`lscol') lwidth(vthin) )
	local words text(`Tht' `lst' "Lean" "Season") text(`cvt' `cov' "COVID-19") text(`Tht' `hvt' "Rice" "Harvest")
	local arrows (pcarrowi `ay1' `ax1' `ay2' `ax2' (12) "2020", color(`cvcol') mlabcolor(`cvcol') mlabsize(medium)) ///
			(pcarrowi `by1' `bx1' `by2' `bx2' (12) "Previous Year", color(gs0) mlabcolor(gs0) mlabsize(medsmall)) ///
			(pcarrowi `cy1' `cx1' `cy2' `cx2' (3) "2019", color(`cvcol') mlabcolor(`cvcol') mlabsize(medium))
	local yaxis ytitle("Index of Household Food Insecurity", axis(1)) ylabel(#5,glcolor(gs12%25))
	local xaxis xlabel(`start'(30.6)`end', angle(45) format(%tdm)) xtitle("")
}

* Make the graph
twoway `shading' ///
	(rarea hi_pre_hunger lo_pre_hunger month, color(gs13) lwidth(none)) ///
	(rcap lo_cvd_hunger hi_cvd_hunger month, lcolor(`cvcol') msize(small) lwidth(medium)) ///
	(line b_pre_hunger month, lcolor(gs0)) (scatter b_cvd_hunger month, mcolor(`cvcol') mlwidth(none))  ///  
	`arrows' , `xaxis' `yaxis' `words' `shade2' ///
	graphregion(color(white)) ylabel(, angle(0))  legend(off) plotregion(margin(zero))
graph export "$dfig/seasonality_graph_`pop'.png", replace		

}
